# 导入必要的库
import numpy as np  # 用于数值计算和数组操作
import matplotlib.pyplot as plt  # 用于绘图和可视化
import seaborn as sns  # 用于绘制更美观的统计图形（如热力图）
from sklearn.datasets import load_iris  # 加载Iris数据集
from sklearn.model_selection import train_test_split  # 用于划分训练集和测试集
from sklearn.cluster import KMeans  # K-Means聚类算法
from sklearn.preprocessing import StandardScaler  # 数据标准化工具
from sklearn.metrics import silhouette_score, adjusted_rand_score  # 聚类评估指标
from sklearn.decomposition import PCA  # 主成分分析，用于降维

# 设置matplotlib的中文字体（解决中文显示问题）
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体字体，支持中文显示
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题，避免出现乱码

# 1. 加载数据集
# 加载Iris数据集，包含150个样本，4个特征（萼片长度、宽度，花瓣长度、宽度）
iris = load_iris()
X = iris.data  # 特征矩阵：150 x 4，包含所有样本的特征数据
y = iris.target  # 真实标签：0, 1, 2（对应setosa、versicolor、virginica），仅用于评估

# 2. 数据预处理
# K-Means对特征尺度敏感（不同特征的量纲可能不同），需要标准化数据
scaler = StandardScaler()  # 初始化标准化器，将特征缩放到均值为0，方差为1
X_scaled = scaler.fit_transform(X)  # 标准化数据，fit_transform同时计算均值和方差并转换数据

# 3. 划分训练集和测试集（可选，用于评估）
# K-Means是无监督学习算法，划分训练集和测试集仅用于评估模型性能
# test_size=0.2表示20%为测试集，random_state=42确保结果可复现
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, test_size=0.2, random_state=42)

# 4. 确定最佳簇数（结合肘部法则和轮廓系数）
inertia = []  # 列表，用于存储每个簇数k对应的惯性（Within-Cluster Sum of Squares）
silhouette_scores = []  # 列表，用于存储每个簇数k对应的轮廓系数
K_range = range(2, 11)  # 定义簇数范围：2到10（k=1对轮廓系数无意义）

# 遍历不同的k值，计算惯性和轮廓系数
for k in K_range:
    # 初始化K-Means模型，n_clusters=k指定簇数，random_state=42确保可复现
    # n_init=10表示运行10次不同初始中心，选择最佳结果
    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
    kmeans.fit(X_scaled)  # 在标准化后的数据上训练K-Means模型
    # 记录惯性（簇内平方和，衡量簇内数据点的分散程度）
    inertia.append(kmeans.inertia_)
    # 记录轮廓系数（衡量簇内紧密度和簇间分离度，范围[-1, 1]）
    silhouette_scores.append(silhouette_score(X_scaled, kmeans.labels_))

# 4.1 绘制肘部法则图
# 肘部法则通过观察惯性随k的变化，寻找拐点（下降速率减缓的位置）
plt.figure(figsize=(10, 6))  # 设置图形大小
plt.plot(K_range, inertia, marker='o')  # 绘制折线图，marker='o'表示数据点用圆点标记
plt.xlabel('簇数 (k)')  # X轴标签
plt.ylabel('惯性 (Inertia)')  # Y轴标签
plt.title('K-Means肘部法则')  # 图标题
plt.grid(True)  # 显示网格线
plt.show()  # 显示图形

# 4.2 绘制轮廓系数图
# 轮廓系数图展示不同k值下的聚类质量，值越大表示聚类效果越好
plt.figure(figsize=(10, 6))
plt.plot(K_range, silhouette_scores, marker='o', color='r')  # 红色折线图
plt.xlabel('簇数 (k)')  # X轴标签
plt.ylabel('轮廓系数 (Silhouette Score)')  # Y轴标签
plt.title('K-Means轮廓系数')  # 图标题
plt.grid(True)  # 显示网格线
plt.show()  # 显示图形

# 4.3 自适应选择最佳k值（基于轮廓系数最大值）
# 轮廓系数最大的k值通常表示最佳聚类效果
optimal_k_silhouette = K_range[np.argmax(silhouette_scores)]  # 找到最大轮廓系数对应的k
print(f"基于轮廓系数的最佳簇数 k: {optimal_k_silhouette}")

# 4.4 自适应选择最佳k值（基于肘部法则的拐点）
# 使用拐点检测方法：计算inertia的二阶差分，拐点通常是二阶差分绝对值最大的位置
inertia_diff1 = np.diff(inertia)  # 计算一阶差分（inertia的变化量）
inertia_diff2 = np.diff(inertia_diff1)  # 计算二阶差分（变化量的变化量）
# 拐点位置：二阶差分绝对值最大的k（+1是因为diff会减少长度）
elbow_k = K_range[1 + np.argmax(np.abs(inertia_diff2))]
print(f"基于肘部法则的最佳簇数 k: {elbow_k}")

# 4.5 综合选择最佳k值
# 优先选择轮廓系数的k值，但如果与肘部法则差异较大（>1），则选择肘部法则的k
optimal_k = optimal_k_silhouette if abs(optimal_k_silhouette - elbow_k) <= 1 else elbow_k
print(f"综合选择的最佳簇数 k: {optimal_k}")

# 5. 训练K-Means模型
# 使用自适应确定的最佳k值训练最终的K-Means模型
kmeans = KMeans(n_clusters=optimal_k, random_state=42, n_init=10)  # 初始化模型
kmeans.fit(X_scaled)  # 训练模型
y_pred = kmeans.predict(X_scaled)  # 获取聚类标签（每个样本所属的簇）

# 6. 模型评估
# 计算轮廓系数，评估聚类质量（范围[-1, 1]，值越大表示簇内紧凑、簇间分离好）
silhouette_avg = silhouette_score(X_scaled, y_pred)
print("轮廓系数 (Silhouette Score):", silhouette_avg)

# 计算调整后的Rand指数，衡量聚类标签与真实标签的相似度（范围[-1, 1]，值越接近1越好）
ari_score = adjusted_rand_score(y, y_pred)
print("调整后的Rand指数 (Adjusted Rand Score):", ari_score)

# 7. 可视化结果
# 使用PCA将4维数据降维到2D，以便可视化（原始4维数据无法直接绘制）
pca = PCA(n_components=2)  # 初始化PCA，降维到2维
X_pca = pca.fit_transform(X_scaled)  # 应用PCA转换

# 绘制聚类结果散点图
plt.figure(figsize=(10, 6))
# 散点图，c=y_pred根据聚类标签着色，cmap='viridis'使用颜色映射，edgecolor='k'添加黑色边框
scatter = plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y_pred, cmap='viridis', edgecolor='k')
plt.xlabel('第一主成分')  # X轴标签
plt.ylabel('第二主成分')  # Y轴标签
plt.title(f'K-Means聚类结果（2D PCA投影，k={optimal_k}）')  # 动态标题显示最佳k值
plt.colorbar(scatter, label='簇标签')  # 添加颜色条，显示簇标签
plt.show()

# 8. 可视化真实标签与聚类标签对比
plt.figure(figsize=(10, 6))
# 散点图，c=y根据真实标签着色，展示真实类别分布
scatter = plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap='viridis', edgecolor='k')
plt.xlabel('第一主成分')  # X轴标签
plt.ylabel('第二主成分')  # Y轴标签
plt.title('真实标签分布（2D PCA投影）')  # 图标题
plt.colorbar(scatter, label='真实类别')  # 添加颜色条，显示真实类别
plt.show()

# 9. 簇中心可视化（降维后）
# 将簇中心降维到2D空间，以便在图中显示
centers_pca = pca.transform(kmeans.cluster_centers_)
plt.figure(figsize=(10, 6))
# 绘制数据点散点图，alpha=0.5设置透明度，显示聚类结果
plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y_pred, cmap='viridis', alpha=0.5, edgecolor='k')
# 绘制簇中心，红色X标记，s=200设置大小，linewidths=3设置粗细
plt.scatter(centers_pca[:, 0], centers_pca[:, 1], c='red', marker='x', s=200, linewidths=3, label='簇中心')
plt.xlabel('第一主成分')  # X轴标签
plt.ylabel('第二主成分')  # Y轴标签
plt.title(f'K-Means簇中心与数据点（2D PCA投影，k={optimal_k}）')  # 动态标题显示最佳k值
plt.legend()  # 显示图例
plt.show()